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^ ! ABSTRACT 

Qh' We investigate the orbital evolution of planetesimals in a self-gravitating circumstellar disc 

in the size regime (~ 1 - 5000 km) where the planetesimals behave approximately as test 
j!^ . particles in the disc's non-axisymmetric potential. We find that the particles respond to the 

' stochastic, regenerative spiral features in the disc by executing large random excursions (up 

Ci ', to a factor of two in radius in ~ 1000 years), although typical random orbital velocities are of 

order one tenth of the Keplerian speed. The limited time frame and small number of planetes- 
imals modeled does not permit us to discern any net direction of planetesimal migration. Our 
chief conclusion is that the high eccentricities (~ 0.1) induced by interaction with spiral fea- 
' tures in the disc is likely to be highly unfavourable to the collisional growth of planetesimals 

\ in this size range while the disc is in the self-gravitating regime. Thus as recently argued by 

Rice et al 2004, 2006, the production of planetesimals gets under way when the disc is in the 
P\| ' self-gravitating regime (either at smaller planetesimal size scales, where gas drag is impor- 

tant, or via gravitational fragmentation of the solid component), then the planetesimals thus 
■ produced would not be able to grow collisionally until the disc ceased to be self-gravitating. 

' It is unclear, however, given the large amplitude excursions undergone by planetesimals in 

, the self-gravitating disc, whether they would be retained in the disc throughout this period, or 

whether they would instead be lost to the central star. 
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1 INTRODUCTION 

It is currently unclear whether the process of planet formation be- 
gins during the earliest phases of star formation (i.e. in circumstel- 
lar discs that are strongly self-gravitating). The high disc mass at 
this stage (with gas mass ~ 10% of the central star's mass (e.g. 
Eisner and Carpenter 2006) is favourable to this hypothesis and 
provides the only opportunity for gas giants to form via (gas phase) 
gravitational instability (Boss 2000). A more subtle effect - again 
strongly favoured during the self-gravitating case - arises through 
the interplay between pressure gradients in spiral arms and gas drag 
on so lid particles and results in the accumulation of dust in spiral 
arms teice et aflliooi) . This dust may then be able to accumulate 
further thr ough the action o f collisions and/or self-gravity in the 
dust phase jRice et alj20o3) .n 

' Note that such rapid accumulation, wfiich might result in the rapid forma- 
tion of large size planetesimals, is required in order for the growing plan- 
etesimals to escape the dangerous meter-size range, where they are sub- 



Once the disc is no longer self-gravitating in the dominant 
(gas) component, it is evidently impossible to form planets through 
either gas phase leans instability or dust accumulation in spi- 
ral arms and thus the most successful mo dels invoke the colli- 
sional growth of dust dPollack et al.l 1 199^ followed possibly by 
the accretion of a gaseous envelope. Despite the lower surface 
densities during later phases (estimates of gas disc mass from 
sub-millimetre dust measurement suggest typical masses that are 
arou nd an order of m agnitude lower than in the self-gravitating 
case; ' Andrews & Wil liams 2005), the great advantage for planet 
formation is simply the longer duration of this later phase. Discs 
with sufficient gas an d dust to form pl a nets typically survive 
for lO'^ - 10^ ye ars faaisch et al] 1200 ll ; lArmitage et d] l2003l : 
IWvatt et all boOS i). whereas the existence of dusty debris discs 
around older stars suggests that (potentially planet building) coUi- 

ject to very fast migration down into the hottest inner disc iWeidenschillind 
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sions between rocky plane tesimals proceeds for 10^ years or more 
jKenvon & Bromlevll2004h . These numbers should be contrasted 
with the brief window of opportunity (~ 10^ years) in the self- 
gravitating phase. 

In summary, then, it is unclear without detailed calculation 
which of these regimes is favoured in the trade-off between lifetime 
and disc mass. 

The viability of planet formation in the self-gravitating regime 
however requires that a number of conditions are satisfied beyond 
the initial phases of accumulation outlined above. In particular, we 
need to understand the dynamical evolution of the protoplanetary 
components (whether at the scale of dust aggregates, planetesimals 
or proto-giant planets) with a view to answering the following ques- 
tions: (i) does orbital migration result in such bodies spiralling into 
the central star, or can they be retained in the disc at least over the 
initial self-gravitating phase, and (ii) (in the case of dust or plan- 
etesimals) are the kinematics of such bodies conducive to further 
coUisional growth? 

As a first step in answering this question we here undertake 
pilot simulations of a small number of test particles in a self- 
gravitating disc, with a view to both considering the issue of their 
orbital migration and to establishing the velocity dispersion that 
is set up through interaction with spiral structure in the disc. It 
should be stressed that the calculation differs in several impor- 
tant ways from most calculations of orbital migration and plan- 
etesimal kinematics reported in the literature. Firstly, conventional 
planetary migration, whether classified as Type I or Type II, is 
the consequence of the gravitational torque exerted on the planet 
(or planetesimal) by non-axisymmetric structure induced by the 
planet/planetesimals. Since this torque contribution can have sig- 
nific ant contributions at the size scale of the Hill radius (or be- 
low: |Bateetai][20o3), it is evidently important that this scale is 
well resolved in numerical calculations, which makes the model- 
ing of the migration of low mass objects particularly challenging, 
especially with codes such as Smo othed Particles Hydrodynam- 
ics (SPH) jde Val-Borro et alll2006h . In the present case, however, 
we are interested in a regime where the disc has pronounced non- 
axisymmetric structure in the absence of the planet/planetesiinal 
(due to the presence of self-gravitating spiral modes in the disc) and 
we are therefore not obliged to resolve the Hill radius in order to 
capture the dominant torque contribution. Evidently, as the mass of 
the planetesimal is increased, the contribution to the torque which 
it experiences from self-induced spiral structure becomes signifi- 
cant and we will assess this a posteriori in order to set an upper 
mass scale on the applicability of our calculations. Secondly, in 
conventional calculations of the kinematics o f planetesimal swarms 
jKokubo &Id"all200(]| : iThommes etZI |2003'), the velocity disper- 
sion of planetesimals is set by an equilibrium between what is 
sometimes termed 'viscous stirring' (i.e. the transfer of kinetic en- 
ergy from larger bodies to the planetesimal swarm via two-body 
scattering) and eccentricity damping by gas drag. In the present 
case, however, the velocity dispersion of our particles is set by the 
level to which eccentricity is pumped through interaction with the 
fluctuating potential of the self-gravitating disc. In this pilot calcu- 
lation we omit gas drag, which places a lower limit to the size scale 
of objects to which the calculation is applicable; again, we set this 
lower limit a posteriori. 

Although we have stressed the qualitative differences between 
this calculation and that reported in the large body of literature on 
planetary migration and planetesimal kinematics, w e point out tha t 
it bears close comparison with the calculation of ( lNelsonll2005h . 
which studies the response of planetesimals to the fluctuating non- 



axisymmetric structure in a disc that is subject to the magneto- 
rotational (MRI) instability. In this study it was found that the plan- 
etesimals underwent stochastic migration, but that it was not possi- 
ble to discern the mean magnitude (or even sign) of the migration 
rate over the limited timeframe of the simulation (a few hundred or- 
bits). In our experiment also, conducted over a similar time frame 
to Nelson (2005), we are unable to discern a net direction of migra- 
tion, but we find (perhaps unsurprisingly given the fact that our disc 
is about an order of magnitude larger in mass) that the amplitude 
of the stochastic variations in the orbital radius is much larger, with 
several planetesimals undergoing radial excursions of a factor two 
or more in both directions. Evidently, this behaviour has important 
implications for the ability of planetesimals to grow by collisions 
during the self-gravitating phase. 

The structure of the paper is as follows. In Section 2 we de- 
scribe the modeling of the self-gravitating disc, which is maintained 
in a 'self-regulated' (i.e. non fragmenting) state through the impo- 
sition of an ad hoc cooling law for which the cooling time is a 
suitably large multiple of the local dynamical time. In Section 3 
we analyse the orbital response of 12 test particles placed in a disc 
of mass 0.1 and O.SMq. We also analyse the location of the domi- 
nant torque contribution on the planetesimals and assess the range 
of planetesimal mass and size scales for which we expect 'test par- 
ticle' behaviour and for which gas drag can be ignored. In section 
4 we discuss the implications of our numerical findings for the mi- 
gration and coUisional growth of planetesimals in self-gravitating 
discs. Section 5 contains our chief conclusions. 



2 NUMERICAL SIMULATIONS 
2.1 Numerical code 

In this paper we simulate the interaction of a gaseous self- 
gravitating disc with a small number of plane t esima l s using SPH , 
a Lagrangian pa rticle based method i ILucvI Il977l : iBenzl Il990l: 
iMonaghM I I992h . This allows the simulati on of gaseous a nd A'- 
body particles, using individual time-steps jBate et al]|l995l) , thus 
resulting in a large saving in computational time when a large dy- 
namic range of timesteps is involved. 

The general setup of our simulations is similar to 
iLodato & Ric"3 ilOOi. l20o3) . We consider a massive gaseous disc 
orbiting around a IMq central star. The disc is taken to extend ini- 
tially from 0.25 to 25 AU and we have considered two cases, where 
the disc mass is equal to 0. IM© (standard runs) and 0.5Mo (massive 
runs), respectively. 

Our simulations employ a number A' = 250, 000 particles and 
are th us intermediate-high resolut ion simulations. Previous simula- 
tions iLodato & Ricel2004l200i) at the same resolution have been 
shown to reproduce reasonably well the internal disc dynamics. In 
particular, once the disc reaches a quasi-steady self -regulated state 
(see below), the disc thickness H is resolved with an average of 2 
smoothing lengths for the low disc mass case and with 4 smoothing 
lengths for the hotter, high disc mass case. 

As mentioned above, our code follows the dynamics of both 
SPH particles and of point masses, which are able to accrete gas 
particles, if they come closer than a given accretion radius to the 
point mass. For the central object, such accretion radius is set to 
0.25 AU. At the end of the simulations typically a fraction no larger 
than 10% of the initial disc mass has been accreted onto the cen- 
tral star. The planetesimals have a nominal mass equal to the de- 
fault gas particle mass, that is 4 x 10"^ for the standard runs 



Eccentricity growth of planetesimals in a self- gravitating protoplanetary disc 3 



and 2 x 10"* Mq for the massive run. This choice ensures that the 
gravitational force of individual planetesimals does not influence 
the overall disc structure and the dynamics of the planetesimals. In 
order to prevent planetesimals from acceting gaseous particles, we 
have set their accretion radius to a very small value. 

The gas particles are allowed to heat up via pdV work and 
artificial viscosity and to cool down with the following simple pre- 
scription 

Au u 
df 

where the cooling time-scale fcooi = /3fl"' and /? is fixed in space 
and time. While this is not a realistic prescription for the cooling 
in actual protostellar discs, it provides a convenient way of param- 
eterising the cooling physics in order that one can ensure that one 
is modeling discs that remain in the self-gravitating regime without 
fragmenting. In par ticular, if the value of the param eter yS is large 
enough (J3 > 5 - 6 dGammidl200lLlRice et"Zl l2005^. it allows the 
realization of a feedback loop where the disc evolves into a quasi- 
steady, self-regulated case, and is maintained close to marginal sta- 
bility. For smaller values of /3 the disc undergoes fragmentation. 
In our simulations we have set /? = 7.5. We note that since this is 
close to the critical value of /J, we expect such a disc to be relatively 
'active' (in the sense that the maintenance of thermal equilibrium 
requires spiral modes to be of relatively large amplitude). We revisit 
this point in Section 4. 




Figure 1. Initial positions of the planetesimals in the evolved quasi-steady 
spiral structure of the 0.1 M© gas disc at the beginning of the run. The 
dimensions are 60 AU in each direction and the logarithmic colour scale 
shows the surface density Z and goes from 10 to 10^ g cm"^. 



2.2 Initial conditions 

Initially the gas in the disc is relatively hot, with a temperature pro- 
file T(R) PC R^^/^ and a minimum value of the stability parame- 
ter ( lToomrelll964) Q - 2, attained at the outer disc edge. Since 
marginal stability corresponds to 2 = 1^ the whole disc is ini- 
tially gravitationally stable. We initially follow the evolution of the 
gaseous disc only, as it cools down and becomes gravitationally 
unstable. The initial evolution is very similar to the one de s cribed 
in several other p apers for example fcodato & Ric3l2004 [20o3 : 
iMeiia et al]|2005l) : when Q ~ I, the disc develops a spiral struc- 
ture, inducing dissipation through compression and mild shocks 
and therefore heating up the disc and establishing a feedback loop 
which keeps the disc close to marginal stability, so that the Toomre 
Q parameter is close to unity over the radial range ~ 1-20. This ini- 
tial evolution takes roughly 60 outer dynamical times for the stan- 
dard runs and 45 outer dynamical times for the massive runs. In this 
state, the surface density distribution is dominated by pronounced 
(of order unity in amplitude) spiral features, which, although regen- 
erative, are individually transient structures (i.e. with lifetimes that 
are of order the local dynamical timescale). 

An ensemble of twelve planetesimals was then introduced be- 
tween 10.8 and 18.8 AU, aligned in four equispaced radial spokes 
(see Figure[T]for the 0.1 M, disc and Figure|2]for the 0.5 M, disc). 
The default velocities were chosen to be Keplerian with a small 
correction for the gravitational potential of the disc, although we 
also investigated the case where the planetesimals were introduced 
with a non-zero initial eccentricity. 



3 RESULTS 

3.1 Orbital evolution in the 0.1 solar mass disc case 

In Fig. [3] we show the time evolution of the orbital radii of all the 
planetesimals for the case where Mrf,,^ = O.IM, and where the 




Figure 2. Same as Fig [T] but for the 0.5 Mq gas disc. Note the warmer 
conditions in the massive disc in the self-regulated state lead to broader 
spiral features than in the standard case (Fig[T) 



planetesimals are initially set in nearly circular orbits. Although 
some planetesimals undergo large amplitude excursions (factor two 
or more), the average orbital radius of the ensemble varies by no 
more than 10% during the simulation, indicating that planetesimals 
in a given region are apparently as likely to go inwards as outwards. 
The thick dashed line in fig. |4] shows the mean eccentricity of the 
planetesimals as a function of time, showing that e ~ 0.05 is typi- 
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Figure 4. Time evolution of the mean eccentricity of the planetesimals for 
various initial eccentricities: e(t = 0) » (bold dashed line), e(t = 0) = 
0.06 (bold solid line), e(t = 0) = 0.03 (thin solid line) e(t = 0) = 0.37 
(thin dashed line). During the course of the simulations a quasi-steady mean 
eccentricity e » 0.07 is established in all cases where the initial eccentricity 
is less than 0.1. 
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Figure 6. Energy and angular momentum (in the inertial frame) are plotted 
for three planetesimals in the massive disc case over the time interval 550 - 
700 yr. The near linear relationship, with the same slope, suggests that these 
three planetesimals are interacting with a mode with pattern speed flp = 




4-5 X IQ- 



, i.e. with co-rotation at ~ 14 A.U.. 



cal (although individual planetesimals may temporarily attain still 
higher eccentricities, ~ 0.3). Thus although much of the fine struc- 
ture in Figure|3]is simply a consequence of elliptical orbital motion, 
it is also evident that in addition some planetesimals undergo abrupt 
and/or large scale changes in orbital radius (for example, one plan- 
etesimal migrates over 20 AU in 4000 years, and another follows a 
period of relative orbital quiescence by migrating ~ 5 AU in about 
200 years). Such migration rates imply radial velocities which are 
at times a significant fraction of the Keplerian speed. 

In Figure |4] we also show the results of starting the planetes- 
imal ensemble with non-zero initial eccentricity, and see that the 
results are unchanged for initial eccentricity less than ~ 0.1. On 
the other hand, if a larger initial eccentricity is employed, the mean 
eccentricity of the ensemble remains around this value for the du- 
ration of the experiment. 



3.2 Orbital evolution in the 0.5 solar mass disc case 

Figure [5] illustrates the rich variety of orbital histories of planetes- 
imals in the massive disc case. The amplitude of orbital migra- 
tion is somewhat larger than in the standard case, although again 
the change in average orbital radius for the ensemble is small, 
< 20%. We note that, as in the standard case, the planetesimals 
undergo stochastic migration events even when they are located at 
radii (> 20 AU) where Q is generally somewhat larger than unity. 
These events coincide with spiral structures temporarily extending 
out into these relatively quiescent regions. 

The mean eccentricity of the ensemble is about double its 
value in the standard case ~ 0.15, with individual planetesimals 
temporarily attaining eccentrcities of up to 0.5. Although some of 
the quasi-periodic structure in Figure[5]is simply attributable to el- 
liptical orbital motion (i.e. at roughly constant angular momentum), 
we also see planetesimals in which, temporarily, the energy and 
angular momentum are subject to correlated quasi-periodic varia- 
tions. Such behaviour suggests that these planetesimals are at this 
stage being driven by a quasi-steady rotating non-axisyrrmietric po- 
tential. We have investigated this possibility by noting that in this 
case the Jacobi constant (Ej) of the planetesimal's motion (i.e the 
total energy measured in the frame of the rotating pattern) should 



be conserved, and that since Ej is related to the energy, E, in the 
inertial frame, and angular momentum, L, via 

Ej = E- LQp (2) 

(where Hp is the pattern speed), we can test for such behaviour by 
plotting E versus L. This exercise indicates that there are indeed 
periods when certain planetesimals appear to be exhibiting this be- 
haviour. This can be evidenced by a linear relation between E and 
L, as shown in Figure (SJ for three of the simulated planetesimals. 
From the analysis of this relation, we can then determine the pattern 
speed of the driving pattern. The planetesimals illustrated in Figure 
[6]appear to be interacting with a mode with co-rotation at around 14 
AU: inspection of an animated series of snapshots of the simulation 
indicates that these three planetesimals are indeed roughly corotat- 
ing with a spiral feature during this period, but that, as the feature 
dissolves shortly thereafter, they then migrate and stall close to an- 
other spiral arm at a different radius. This behaviour is not seen in 
the standard (0.1 Mq disc) simul ations, which is cons istent with the 
expectation (see, for example, Lodat o & Ricell20 0^ that the more 
massive disc produces rather more coherent and long-lived spiral 
features. Nevertheless, we find only occasional evidence for plan- 
etesimals that are being driven by a dominant quasi-steady spiral 
mo de so that the analysis of orbita l families in this situation (see 
e.g. lVorobvov & Shcheki"novll2006h is of limited applicability here. 

3.3 Planetesimal-planetesimal encounters and gas densities 
close to the planetesimals 

Inspection of animations of the simulations suggest that planetes- 
imals tend to spend time preferentially in spiral structures in the 
disc. As such structures dissolve on a roughly dynamical timescale, 
the planetesimals then migrate and find further temporary lodging 
in a new spiral structure. (The strongly fluctuating nature of the 
spiral potential means that, as discussed above, the planetesimals 
usually have insufficient time to attain steady orbits in the frame 
co-rotating with the local spiral pattern). This predilection for be- 
ing in the spiral arms suggests that the planetesimals are expected 
to come closer to each other and to sample higher gas densities than 
in the case of particles orbiting in an axisymmetric disc. 

We illustrate the increased tendency for planetesimal- 
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Figure 3. Evolution of the distance to tlie star of three representative planetesimals for the M^i^ = OAMq case. While individual planetesimals might migrate 
significantly over the course of the simulation, the average orbital radius only varies by less than 10%. 
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Figure 5. Evolution of the distance to the star of three representative planetesimals for the Mdisc = O.SMq case. 
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Table 1. Number of encounters with closest approach smaller than d (mea- 
sured in AU) between planetesimals during the first 15000 time units in the 
0. 1 M» disc run and the non-self-gravitating disc run. 



planetesimal encounters in Table [T] where we contrast the num- 
ber of encounters between the planetesimals in various separation 
ranges with the corresponding statistics from a control run in which 
the planetesimals orbit in a smooth, non-self-gravitating disc. We 
see that the number of encounters with closest approach smaller 
than 0.5 AU is significantly enhanced in the self-gravitating disc 
case. As expected, the effect of planetesimals lingering near spiral 
features is to decreasse the minimum distance between planetesi- 
mals in the simulation. We stress that, given the very low masses 
of the planetesimals (i.e. equal in mass to a single gas particle), the 
effect of mutual encounters on the planetary dynamics is negligible 
- as intended - even for the closest encounters in the simulation. 

We illustrate the effect of non-axisymmetric structure in the 
disc on the gas density field sampled by the orbiting planetesimals 
in Fig. |7] We plot - for a given planetesimal - the density excess, 
which we define as 



2(Pp -Pr) 
Pp+Pr 



(3) 



where pp is the gas density in the vicinity of the planetesimal and 
Pr is the azimuthally averaged density at the instantaneous orbital 
radius of the planetesimal. This planetesimal appears to spend most 
of its time in the spiral arms and thus at larger gas densities. This 
is typical, although some planetesimals also undergo phases where 
they sample lower gas densities as can be seen in Fig. ??. Such 
periods in under-dense regions are less apparent in the case of the 
massive disc simulation. 

We see from Figure 7 that although this planetesimal samples 
gas densities larger than the local azimuthal average, 77 is always 
less than or of order unity. This is to be expected: the amplitude of 
the spiral features in the gas is only of order unity and thus even if 
planetesimals spent all their time in spiral arms, they would only 
experience a density that exceeded the local mean by a modest fac- 
tor. 

3.4 Numerical issues 

Firstly, our experiments have aimed to study the response of each 
planetesimal to the fluctuating disc potential and have only mod- 
eled the evolution of an ensemble of planetesimals for reasons of 
computational economy. Therefore, we need to be satisfied that the 
stochastic behaviour we see is a result of planetesimal-disc and 
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Figure 7. Comparison between the density in the vicinity of one planetes- 
imals in the 0. 1 M, disc with the mean gas density at the current radius of 
the planetesimal. 
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Figure 8. Evolution of the density excess (see equation (3) as in Figure 
7, but for a different planetesimal. 



not planetesimal-planetesimal interactions. We find, in fact, that 
the torques experienced by the planetesimals due to other plan- 
etesimals are typically two to three orders of magnitude less than 
those arising from the disc. This is to be expected. The closest 
planetesimal-planetesimal interactions occur at distances of ~ 1 
AU: given that a typical smoothing length in the disc is ~ 0.2 AU 
and that a smoothing kernel contains 50 particles, it follows that 
even at closest approach there are about 6000 gas particles that are 
closer to a planetesimal than its nearest planetesimal neighbour. 
Given that the planetesimal masses in the simulations are equal to 
the gas particle masses, it follows that the gravitational torques ex- 
erted between planetesimals are negligible. 

Secondly, we are interested here in examining the regime 
where the dominant torques experienced by the planetesimal de- 
rive from the non-axisymmetric density distribution that the disc 
would have in the absence of the planetesimal (i.e. that due to the 
spiral structure in the disk) rather than the non-axisymmetric struc- 
ture induced in the disk by the planetesimal. Note that in non-self 
gravitating discs, and thus in all discussions of planetesimal migra- 
tion in the literature (with the exception of that in Nelson 2005) it is 
the latter effect that is responsible for planet migration. In this case, 
it is essential that the calculation properly resolves the Hill radius 
of the planetesimal 
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Figure 9. Integrated z component of the torque as a function of distance 
from the planetesimal, normalised to the torque from within 20 AU of the 
planetesimal. 



where is the radial coordinate of the planetesimal and Mp the 
mass of the planet, since this is the region from which the bulk of 
the torques on the planetesimal arise. 

In the present case, by contrast, where we examine migration 
driven by spiral structure in the self-gravitating disc, it is no longer 
important to resolve the Hill radius and we can thus examine the 
migration of low mass planetesimals for which resolution of the 
Hill radius would be impracticable. Nevertheless, we need to check 
that our assumption is correct, i.e. that the dominant torques indeed 
arise at relatively large distances from the planet. Fig.|9]illustrates a 
typical snapshot of the cumulative distribution for the z component 
of the torque of the disc on the planetesimal as a function of dis- 
tance from the planetesimal. We see that the dominant contribution 
to the torque is around 1 AU from the planetesimal and that the con- 
tribution to the torque from the poorly resolved region around the 
planetesimal (i.e. from within a typical particle smoothing length in 
that region, around 0.2 AU) is small (< 20%). We are thus satisfied 
that the orbital evolution is not driven by the disc's response to the 
planetesimal and that the planetesimal should be behaving approx- 
imately as a test particle. We examine this assumption further by 
investigating how the migration depends on planetesimal mass in 
Section[331 

On the other hand, in the case of the component of the torque 
perpendicular to z, i.e. acting on the inclination, the contributions 
from close to the planetesimal actually can make up a large frac- 
tion of the torque (i.e. around 40% of the torque originating from 
within a smoothing length). Thus the inclination development of 
the simulations cannot be trusted. Nevertheless the inclinations stay 
relatively small (typically less than 0.02 radians in the standard and 
0.03 for the massive disc). This is small compared with the axis ra- 
tio of the disc H/R ~ O.I and so the planetesimals remain mainly 
confined within the disc. 

3.5 Valid regime 

We now discuss the range of planetesimal masses for which the cur- 
rent calculation is valid. As noted in Section [2n we are modelling 
the regime in which the planetesimal behaves as a test particle in 
the gravitational potential of the star-disc system and where the re- 
sults should be independent of planetesimal mass. In reality, at low 
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masses the effect of gas drag (omitted in the present simulations) 
would break this mass independence, whereas at high masses, the 
disc response to the presence of the planetesimal affects the torques 
on the planetesimal. We consider each of these issues in turn. 

A lower limit for the mass range of our simulation can be ob- 
tained by calculating the mass scale at which the stopping time due 
to drag is comparable with the timescale on which the planetesimal 
undergoes stochastic migration due to interaction with spiral struc- 
ture in the disc. From inspection of Figure[3]we see that planetesi- 
mals typically reverse their direction of migration on timescales of 
about 10' years. Thus if we demand that the timescale for gas drag 
to significantly perturb the orbit should exceed say 10* years (~ 58 
orbits), we should safely be in the regime where gas drag plays a 
minor role in the dynamics of the system. 

The drag force on solid bodies within a pro tostellar disc Fo is 
given by jWeidenschillindl 19771 : Iwhippldl 1972h 

= ^CD7raV"^ (5) 

where a is the radius of the planetesimal, u is the velocity of the 
planetesimal relative to the gas, p is the gas density and Cd is the 
drag coefficient. 

Thus the timescale on which relative motion between the plan- 
etesimal and the disc gas is damped by gas drag is 



(6) 



In the case that we are considering here, the relative velocity be- 
tween the planetesimal and the disc gas does not derive - as in the 
case usually considered (e.g. Weidenschilling 1977) - from the dif- 
ference between Keplerian planetesimal motion and sub-Keplerian 
gas motion subject to outward pressure forces. Instead, this relative 
velocity, u, derives from the eccentricity of the planetesimal orbits 
induced by the spiral structure in the disc, so that we have u = eV^, 
where Vk is the local Keplerian velocity. 

Assuming a spherical planetesimal of density pp, we can then 



write t. as 



fe = 



SppR 



(7) 



or, equivalently (and adopting Co in the range 0.5 - 1 , which is 
appropriate to the case considered here where the Reynolds number 
exceeds a few hundred) 



(8) 



8 Pp a Q-^ 

where R is the orbital radius of the planetesimal. If we then require 
that fe > 10"* years (~ SOOfl"'), and adopting typical values of 
e ~ 0.1, p ~ a few xlO"" g cm"' and pp ~ a few g cm"', we then 
obtain the requirement that a must be at least of order a kilometre 
or so before one can neglect the effect of gas drag. 

We thus conclude that the effect of gas drag is negligible for 
planetesimals of kilometre scale or above (corresponding to bodies 
of masses ~ 10'* g.) 

Our simulations are valid for all higher masses as long as the 
accelerations experienced by the planetesimals remain independent 
of mass. But these accelerations change as soon as the planetesi- 
mals have enough mass to gravitationally influence the gas density. 
Thus we get an upper limit for the validity of our simulations by 
checking at what mass scale the torque per planetesimal mass T/Mp 
changes with mass. In Figure 10 we plot the evolution of T/Mp for 
a single planetesimal in the case that it has the mass indicated in 
each of the curves. It can be seen in Fig.[lO]that significant changes 




Figure 10. The torque per unit planetesimal mass as a function of time for 
runs in which the mass of the planetesimal is respectivefy 1 / f (bold dashed 
line), 1 (bold solid line), 10 (thin solid line) and 100 (thin dashed line) times 
the mass of a gas particle. 
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Figure 11. The distance to the star is plotted against time for one planetes- 
imal in the normal run (solid line), the run with a tenth planetesimal mass 
(long dashed line), the run with 10 times higher planetesimal masses (short 
dashed line) and the run with 100 times higher planetesimal masses (very 
short dashed line). 



to the torques on the planetesimals only arises when the planetesi- 
mal mass is raised above ten times the one used in our standard run, 
corresponding to about 8 x 10^' g. We reach the same conclusion 
from inspection of Fig.[TT] where again we see a marked deviation 
in orbital evolution for the largest mass planetesimal. 

So summarising, we have as a region of validity of our simu- 
lations: 



2 • 10"** M„ 



10 



" < M < 8 ■ 10^' g = IMe 



(9) 



(10) 



or in planetesimal radii 
Ikm <a< 5000km. 

4 DISCUSSION 



We can summarise our numerical findings as follows. Planetesimals 
in the size/mass range to which our calculations apply (see equa- 
tionllOll are subject to large scale fluctuations in orbital radius due 
to interaction with regenerative spiral structure in the disc. We can 
discern no net sign to the mean migration rate over the timescale 
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of the simulation and are thus unable to answer one of the goals of 
the present study, i.e. whether we expect planetesimals to migrate 
into the star during the self-gravitating phase of the disc. Since the 
disc's lifetime in the self-gravitating phase exceeds the duration of 
the present experiment by about a factor 30 it will be computation- 
ally challenging to answer this question, especially since (given the 
stochastic nature of the orbital evolution) it will be necessary to 
assess the fraction of planetesimals retained in the disc through fol- 
lowing the evolution of a large ensemble of particles. 

We can however make some comments about how the dy- 
namical evolution of our particles is likely to affect their col- 
lisonal growth during the self-gravitating phase. As noted by Nel- 
son (2005) in the context of the evolution of planetesimals in a disc 
subject to MRI, the fluctuations in orbital radius and eccentricity 
could in principle either aid or hinder planetesimal growth. On the 
positive side, orbital migration provides a mechanism for avoid- 
ing the self-limitation of planetesimal growth through the clear- 
ing out of a proto-planet's "feeding zone". In conventional mod- 
els, this imposes an "isolation mass" limit which makes it difficult 
to gr ow b eyond a terrestrial mass scale in the inner regions of the 
disc jlda & Lin 2004). Although this limitation is by-passed in our 
self-gravitating models, since the feeding zone would be continu- 
ously replenished by planetesimals undergoing an orbital random 
walk, this factor is not very significant since in any case the high 
surface density during the self-gravitating phase means that isola- 
tion masses are high (higher, that is, than the upper limit for the 
applicability of our calculations). Another potentially positive as- 
pect of evolution in a self-gravitating disc is the possibility of con- 
centrating planetesimals in high density spiral features, where they 
might be more prone to undergo collisional growth. In the case of 
the more massive - O.SMq - disc (see Section 4), there is some 
evidence for the dominance at times of spiral modes which persist 
over a number of dynamical times, allowing particles to settle into 
orbits which conserve a Jacobi integral in the co-rotating frame. 
Some of these orbits involve a preference for particles to spend 
time in the spiral arms, an impression strengthened by inspection 
of animations which show particles which co-rotate for a while in 
the vicinity of an arm and which then respond to changes in the 
spiral structure by passing rapidly to another spiral feature. This 
tendency will not lead to a significantly enhanced collision rate, 
however, unless there is a strong enhancement in the solid to gas 
ratio in the spiral arms: since the amplitude of the gaseous arms is 
around a factor two at most, then even a particle that spent all its 
time in the arms would only, at fixed solid to dust ratio, experience 
an enhancement in the local planetesimal density of a factor two 
at most. Since we have only studied the evolution of ten planetesi- 
mals, we cannot, on the basis of the present study, comment on the 
degree to which planetesimals ar e differen tially concentrated in the 
arms. Nevertheless, the results of lRice etal. (200 4) suggest that the 
diiferential accumulation of solids in the arms is weak at size scales 
where gas drag is unimportant. 

We now turn to the negative implication of stochastic migra- 
tion noted by Nelson (2005), i.e. the growth of planetesimal veloc- 
ity dispersion. This hinders planetesimal growth in two respects: 
higher collisional velocities are associated with disruption of parent 
bodies instead of collisional growth (e.g. Leinhardt et al 2007) and 
also suppress the role of gravitational focusing in effecting physi- 
cal collisions. Quantifying the second effect, the timescale for col- 
lisional growth of a planetesimal of mass Mp, density p and radius 
Rj, is given by: 



^'"^ Ipa(I-H4GMp/(iJpOr2)) 

where Ep is the surface density of planetesimals, cr is the planetes- 
imal velocity dispersion and Q. is the local Keplerian orbital fre- 
quency. For a planetesimal of density ~ 3 g/cm' located at 10 A.U. 
in a disc with surface density of planetesimals equal to 0.4 g/cm^, 
the growth timescale via two body collisions is: 

?g„„ = 5 X 10'- , yrs (12) 

1 -I- M^' V"- 

where M = (Mp/6 10"" Me) and a = (cr/lkm s"'). Evidently this 
timescale is far too long to be of interest unless the second term on 
the denominator (due to enhancement of the collisional cross sec- 
tion by gravitational focusing) is large. If we simply set cr to be the 
product of the orbital eccentricity and local Keplerian velocity (an 
assumption we revisit below), it follows that at a given location in 
the disc and for a given mass scale of planetesimal, then the growth 
timescale in the gravitationally focused regime scales as the square 
of the eccentricity divided by the surface density normalisation. In 
other words, for a planetesimal at a given mass scale to be able to 
grow significantly over a disc phase of duration t\[f^ we require 

^ < 1 (13) 

?life 

which translates into a lower limit on the product t^fj^e"- (where X 
is the surface density normalisation of the disc). 

We can now return to the issue raised in the Introduction of 
whether or not planetesimal growth is favoured in the high density 
but short lived phase during which the disc is self-gravitating. As- 
suming for simplicity that the self-gravitating phase lasts ~ 10% 
of the subsequent lifetime of the primordial gas/dust disc but that 
its surface density (mass) is typically ten times greater, we see 
that these two factors roughly cancel. Thus whether the growth 
of planetesimals is favoured during the self-gravitating phase ac- 
tually boils down to the relative values of the eccentricity in the 
two phases. 

Our simulations have shown that this factor is decisive in 
disfavouring the self-gravitating regime. In conventional mod- 
els for the dynamical evolution of planetesimal swarms (e.g. 
Kokubo & Ida 2000) the equilibrium eccentricity (attained through 
the balance between stirring by larger bodies and damping by gas 
drag) is very low (of order 0.01). We have however seen that in the 
self-gravitating phase, the typical eccentricity is 0. 1 or more. Since 
the growth time scales as the square of the eccentricity, this means 
that the self-gravitating phase is highly unfavourable (compared 
with the later more quiescent phase) for the collisional growth of 
planetesimals. We note that the above estimates for conventional 
proto-planetary discs assume a laminar disc structure; eccentricity 
driving by interaction with turbulent structures in a disc exhibiting 
MRI would further hinder collisional growth, whether or not the 
disc is self-gravitating. 

There however remain two possible qualifiers of the above 
conclusion. Firstly, it is only valid to set the local velocity disper- 
sion equal to the product of the eccentricity and local Keplerian 
velocity in the case that the kinematics is well described as a set 
of elliptical orbits which are elongated in a random direction. Such 
a description does not of course allow for the possibility of local 
velocity coherence in a planetesimal swarm in which velocities are 
significantly non-circular. In this sense, therefore, the above esti- 
mates are pessimistic. Further calculations, involving a large swarm 
of planetesimals from which one could derive the local velocity 
dispersion, are required in order to settle this issue. We however 
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note that our difficulty in identifying coherent orbital families, in- 
teracting with a long lived spiral mode, probably implies that this 
effect is unlikely to improve the prospects for particle growth sig- 
nificantly. The second possible qualifier is that, as noted in Section 
2, our simulations employ a cooling timescale which is close to the 
critical one for disc fragmentation and therefore represent condi- 
tions where spiral modes are of relatively high amplitude. In the 
inner regions of protoplanetary discs, however, cooling timescales 
are longer than employed here and thus one would expect lower 
amplitude spiral structures, with a correspondingly lower ampli- 
tude of stochastic driving of planetesimal orbits. This possibility 
needs to be quantified through further numerical investigation. 

If one leaves aside these possibilities, then we have shown that 
the self-gravitating phase is unfavourable for collisional growth in 
the regime of planetesimal scale studied. We have not ruled out 
the possibility of collisional growth at smaller scales (where gas 
drag can effect a signi ficant enhanceme nt of the local solid to gas 
ratio in spiral features: iRice et al]|2004b nor that (where such en- 
hancement exists) it may not be possible to produce much larger 
bodies throu gh the gravitatio nal fragmentation of the solid com- 
ponent (e.g., i Rice et ai] |2006h . We nevertheless consider it likely 
that whatever the size scale attained in this way, there will be no 
further collisional growth of planetesimals while the disc is self- 
gravitating. 

This obviously limits the extent to which planet formation can 
get under way in the self-gravitating phase. Either a massive planet 
is produced via gas phase Jeans instability CBoss..2000.) , thus obvi- 
ating the need for collisional growth , or else, if self-gravity aids the 
initial accumulation of solid bodies jRice et alj2004l |2006|) . then it 
is unlikely that such bodies can undergo further collisional growth 
once they become large enough for gas drag to be unimportant. In 
this latter case, further growth of such bodies could only resume 
once the disc had evolved to the point where its self-gravity (and 
the associated pumping up of the planetesimal velocity dispersion) 
has become unimportant. The vital unanswered question, therefore, 
is whether orbital migration would prevent the planetesimals from 
remaining in the disc throughout the self-gravitating phase. This 
question can only be answered through further, computationally 
challenging, calculations. 
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